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MTHSGSOQRT ; G Floating Point Square Root routine 16-SEP=1 AX/VMS Macro V04-00 P 1 
1-004 ° 7 g-8 EP=-19 382 9 ‘31 :36 EMTHRTL. SRE IMTMGSORT MARS 1 a8 (1) 


-TITLE MTHS$GSOQRT 3 G fioeting Point Square Root routine 
IDENT /1-004/ : File: MTHGSQRT.MAR EDIT: RNH1004 


eee 


COPYRIGHT (c) 1978, 1980, 1982, 1984 B 
DIGITAL EQUIPMENT CORPORATION, MAYNARD. MASSACHUSETTS. 
ALL RIGHTS RESERVED. 


THIS SOFTWARE IS St | UNDER A LICENSE AND MAY BE USED AND COPIED 
ONLY IN ACCORDANCE WITH THE TERA S OF SUCH LICENSE AND WITH THE 
INCLUSION OF THE ABOVE pd NOTICE. ad SOFTWARE OR ANY OTHER 


THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT NOTICE 
cORPORAT ion NOT BE CONSTRUED AS A COMMITMENT BY DIGITAL EQUIPMENT 


DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY OF ITS 
SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL. 
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of ACeLltys MATH LIBRARY 
"ABSTRACT: 


MTHS$GSQRT is a function which returns the G-floating square root of 
its G-floating argument. The call is standard call-by-reference. 


MTH$SGSQRT_RS is a special routine which is the same as MTH$GSQRT 
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0 except a faster non-standard JSB call is used with the argument in 

3 RO and no registers are saved. 
00 -- 
00 | 
00 VERSION: 1 
00 
HISTORY: 

AUTHOR: ‘ 

bp Steven B. Lionel, 08-Jan-79: Version 1 

0 MODIFIED BY: 

00 

00 

00 


Y ; Detailed Current 
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MTHRTL. Bae SATH HGSQRT.MAR; 1 
History 


\ 
ated. Results should be 


SBL 08-Jan-79 


t 
Use ASHQ, not ee to generate reserved operand. Je Smee 


lo 
R 
; -SBTTL HISTORY ; Detailed Current Edit 
5 
§ 3; ALGORITHMIC DIFFERENCES FROM FP=11/C ROUTINE: 
3 wort used for comparison, FP-11 has no G-float. 
3 : . Last DIVD is rounded instead of trunc 
9 26 correct within 2 LSB's 
+44 °) : Edit History for Version 1 of MTHS$GSQRT 
0000 6 : 1-001 = Adapted from sense version 1-010. 
44 ¢: : be - Correct some os in comments. JBS 30-JUL-1979 
0000 66 : 1-004 - Change w* to G* on call to MTH$SSIGNAL RNH 09-Sept- 1 


MTHSGSQRT ; G Floating Point Square Root ro mo 16-SEP-1984 01: AX/VMS Macro Vv04-00 P 
1-004 BECLARATIONG 3 Declarative Part oy Modul - nets 9 $3! 3§ YRTARTL SRC JMTHGSOQRT.MAR; 1 — (3) 
$3 -SBTTL DECLARATIONS 3; Declarative Part of Module 
o 3 
4 ; INCLUDE FILES: OTSPARAMS .MAR 
0 rs > EXTERNAL SYMBOLS: 
00 75° .DSABL GBL : Declare all externals 
00 4 ~EXTRN mTHSSS SIGNA 3; SIGNAL * SEVERE error 
00 7 -EXTRN MTHSK » SOUROONEG Error code 
8 
$8 80 : EQUATED SYMBOLS: 
000 81; 
$000 8S 
0000403C 444 az ACMASK = “M<IV, R2, R3, R4, RS> ; register save mask and IV enable 
0000 85 > MACROS: none 
0000 6 ; 
44 He ; PSECT DECLARATIONS: 
00000000 89 -PSECT _MTHSCODE PIC,SHR,LONG, EXE ,NOWRT 
0000 90 ; program section for math routines 
0000 91 ; 
0000 3 > OWN STORAGE: 
0000 93 ; 
0000 9% ; 
0000 95 ; Constants A and B chosen for k = odd 
0000 96 ; 
0000 97 LG_ODD_A -f511: 
0000 A000 8279 SFFA 0000 98 WORD ‘%0057772, *0101171, *0120000, *0000000 
0008 99 LG_ODD_B_ “EMS11: 
0000 4000 0789 2903 is 199 ~WORD *0020003, *0003611, *0040000, *0000000 
0010 198 : Constants A and B chosen for k = even 
0010 10 
0010 104 iG. EVEN_A: 
0000 4000 BEC3 4002 at 105 WORD %*0040002, *0137303, *0040000, *0000000 
0018 106 LG_EVEN_B _EMS12: 
0000 6000 £964 1FFA 0018 #107 -@ORD *9017772, *0164544, *0060000, *0000000 


MTH' 
1-0 


a + el int Square Root rout ine 1 mite ts ve 91: 3)! 36 yaalyas Be Macro V04-00 


andard G-Floating GSQRT “SEP=19 
~SBTTL MTHSGSQRT = Standard G-Floating GSORT 
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{FUNCTIONAL DESCRIPTION: 

; GSQRT = G-floating point function 

: GSQRT(X) is computed using the following approximation technique: 
If X <0, error. If X = 0, return GSQRT(X) = 0. 

Let X = 2**K * F where F is the fractional part. 


= Qee( 


2P) 
CORT OX) = 2eeP « *GSORTCFD, 1/2 =< F <1 


If K = odd, X = 2ee(2Pe1) * F = 2ee(2P+2) » (F/2) 
SORT (x) 2e2(P+1) ® GSORT(F/2), 1/4 =< F/2 < 1/2. 


Let F* = AeF + B, 


If K = even 


A = 0.453730314(octal), 
B = 0.327226214(octal), for K = even. 
= At(F/2) + B, 
A= 8 go 8T er egtoctal). 
4 B = 0.230170444(octal), for K = odd. 
an 
K' =P, for K = even 
=P +1 for K = odd. 


Let YO = 2**k' * F* as a straight Line egorey inet ion within the 
given interval uate coefficients A and B which minimize the 
absolute error at the midpoint and endpoint. 


YCn+1) = (1/2) * ( Yn) + X/¥CEn]) 
; The relative error is < 10**-17. 
; CALLING SEQUENCE: 

; gsqrt.wg.v = MTHS$GSORT(x.rg.r) 
INPUT PARAMETERS: 


LONG = 4 3; define Lenquory multiplier 
x = 1 * LONG ; Contents of x is the argument 


IMPLICIT INPUTS: none 
OUTPUT PARAMETERS: 

VALUE: G-floating square root of the argument 
IMPLICIT OUTPUTS: none 
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MTHRTL.SRC JRTHGSORT.MAR; 1 


Starting with YO, three Newton-Raphson iterations are performed. 


Bl, 
int Square Root routine 1 mi tabs 9 m1: 3) 36 me F o V04-00 Page 5 
andarad G-Floating GSQRT -SEP-1984 4:2 MTHRTL.S CIMTHGSORT MAR: 1 (4). 


COMPLETION CODES: none | 
SIDE EFFECTS: | 
Signals: MTHS_SQUROONEG if xX < 0.0 with reserved porate in RO/ 


(copied to thé signal mechanism vector CHFSL_MCH_RO/R1 hy CIBSSTGNAL). 
Associated message is: “SQUARE ROOT OF NEGATIVE ALUE"'. Re 


MTHSGSOQRT ; _G 
1-004 H 
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sult is reserved | 
operand -0.0 unless a user supplied (or any) error handler changes CHFSL _MCH_RO/R1 


NOTE: This procedure disables floating point underflow, enables integer 
overflow, causes no floating overflow or other arithmetic traps, and 
preserves enables across the call. 
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sta 
; disable DV (and FU), enable IV 
MTHSFLAG_ JACKET ; fla 


MOVAB G*MTHSSJACKET_HND, (FP) 
; set handler address to iacket 
3; handler 


Q that this is a jacket procedure in 
6D O0000000'GF 9E 


18 3 case of an error in special routine 
186 MOVG ax(AP), RO ; RO/R1 = ar 

187 BSBB MTHSGSQRT_RS ; call kernel GSQRT rountine 

188 T 3 return - result in RO/R1 
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«ENTRY MTHSGSQRT, ACMASK 3 ndard call-by-reference entry 
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MTHSGSORT ;_G Floating Point Square Root routine 16-SEP-1984 01:31:5 AX/VMS Macro v04-00 P ATH! 
1-004 ATASGSORT RS 2 Spectat GsaRt routines  ‘SoSEPa}obe 44:24:36 EMTMRTe eaeaaTNesans mars? 28° (8) 1-0 
! 139 -SBTTL MTHSGSQRT_RS = Special GSORT routine 
} 136 : Special GSQRT - used by the standard routine, and directly. 
0 1 194 ; CALLING SEQUENCE: 
1 195; save anything in R2:R5 
031 196 : MOVG ase ; input in RO/R1 
031 197: JSB MTH$GSORT_RS 
3 1 198 ; return with resulf® in RO/R1 
1 199 ; Note: This routine is written to avoid causing any integer overflows, floating 
0 : a} 3 over tious, floating underflows or divide by 0 conditions, whether enabled or | 
3 not. 
1 3 
03} 508 REGISTERS USED: | 
031 04 ; RO/R1 = Floating argument then result 
0031 05 ; R2/R3 = scratc ; 
83 1 36 ; R4/R5 = hold X during calc of F", K*. 
0031 208 MTHSGSORT RS: ; JSB routine for GSQRT | 
54 50 SOFD 0031 09 MOVG RO, R4 3 test sign of X and save it in R4/R5. 
54° «O15 0033 10 BLEQ  ZERO_NEG : branch to ZERO_NEG if X =< 0 
Bi Sei > 0 
Se ee $833 14 Pos: MOVZWL RO, R2 ; isolate low 16 bits in R2 
OO3A 15 i; (sign, exp, 4 fract bits) 
38 1F BA QO3A 18 BICB2 #*X1F, R2 ; clear fraction, 1 exp bit 
5 52 AA 003D 1 BICW R2, R ; clear exponent bits in RO 
0c 50 04 €1 sh 18 BBC #4, RO, EVEN ; branch if exponent even 
50 B88 AF 44FD 0044 305 MULG2 LG.ODD.A_E511, RO ; add 511 half of bias) to 
; (exponent- 
0049 1 3 ond’ start approximation calc | 
50 BB AF 40FD Bee ¢ ADDG2 LG_ODD_B_EM511, RO ; RO = (first approx) * 2**-512 
OF 11 4 BRB ADJUST 3 go adjust 
O03) Se ey 
50, 2000 BF Ad 0080 226 “ appw2 -#*x2000, RO ; exp is 0 - make it 512 | 
B7 AF 44FD 05 MULG LG_EVEN_A, RO 
BA AF OFD 035A 228 ADDG2 © LGEVEN B_EMDI2, RO ; RO = (first approx) * 2#*-512 | 
52. 52 «1F «69C «(005F 0 " ROTL #31, R2, R2 ; divide R2 (exptbias) by 2, | 
et: 1 3 giving (exp/2+512) 
50 52 a0 006 § ADDW R2, RO ; insert exp/2 in first approx and 
898 ? 3 re-bias it. | 
04 5 ; first iteration, all done in G floating because of exponent range 
° | 
52 54 50 47FD 066 $ : DIVG3 RO, R4, R2 ; R2 = | 
50 3 40FD 0068 8 ADDG2 R2, RO ; RO = YO + x/Y0 
a Ae ‘sy 9 SUBW #*x10, RO ; RO = ¥1_= (1/2)(YO + X/Y¥O) 
fi rt 3; no overflow possible | 
34 45 3 second iteration 
52. 54 50 47FD 007 44° DIVG3 RO, R4, R2 ; R3/RS = x/Y1 
50 5 40FD 7 45 ADODG Re, ; RO/R1 = Y1 + higher pecs (acy) 
50 1 A2 78 46 SUBW #*x10, RO 3; RO/R1 = Y2 = (172) TY1+x/¥1) 
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MTHSGSQRT Floatin -S 1:31:5 AX/VMS Macro V04-00 
1-004 rise GSORT_RS = Special ESRT routine bE 1984 iE 34 35 EATHRTL. SRC JMTHGSORT.MAR; 1 (5) | 
7 47 | 
f rt ; third iteration 
52. 54 5 47FD 007 . DIVG Rg, Rs R4, R2 : R2/R3 = K/V3 
8 — ADDG ; RO = Y24x/Y 
1 A SUBW Be e16° RO ; RO/R1 = GSQRT(X) = 
A 3 (172) (¥24xK/Y¥2) 
05 : : SQRTX: RSB 3: return, RO/R1 = result 
§ 7; x =< 0 
8 . 
B : 2ERO_NEG 
FD = s«13 8 SQRTX 3 return with RO = result = 0 
ge ge D 60 PUSHL (SP) 3 return PC even JSB routine 
7 00' 8F A F 61 MOVZBL #MTHSK_SQUROONEG, -(SP) ; ROvAt = on value 
50 1 79 : ; $6 ASHQ #15, #T, RO 3 4 /R1 = result = reserved operand -0.0 
9 6 ; RO/R1 goes to signal mechanism vector 
097 64 3 Pearse MCH_RO/ Ri so error handler 
097 65 3 can modify the result. 
00000000'GF 02 097 66 CALLS #2, G*MTHSS$SIGNAL 3; signal error and use real user's PC 
009E 67 3 independent of CALL vs JS 
4943 o8 RSB ; return - RO restored from CHF$L_MCH_RO/R1 
009F % ~END 
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MTHSGSQRT i os Sat 4 9 91:3): 35 yaa /yes Oe o V04-00 Page 
Symbol table 6-SEP-1984 11: MTHRTL.S SMTHGSORT.MAR; 1 
ACMASK = 403C 
ADJUST 1 R 1 
EVEN R 1 
LG_EVEN_A 10 R 1 
LG_EVEN B Ente 18 R 1 
LG-ODD_A_E511 0 K 1 
L& 0008 ER511 0 ; R 1 
MTHSSJACKET_HND eeeeeere x 1 
MTHSS$SIG eeeerere § 0 
MTHSGSQRT 00000020 RG 1 
MTHSGSORT_RS 00000031 RG 1 
THSK_S eeeerene x 0 
POS 80098 37 R 01 
SQRTX OO000BA R 01 
x = 00000004 
ZERO_NEG 00000088 R 01 
torn me cmon se eee seo 
: ; Psect synopsis : 
PSECT name Allocation PSECT No. Attributes 
S 00000000 ( 0.) 00 ¢ O.) NOPIC USR CON ABS LCL NOSHR NOEXE NORD NOWRT NOVEC BYTE 
MTHSCODE OOOOOO9F (¢ 159.) O01 ¢ 1.) PIC USR CON REL LCL SHR EXE RD NOWRT NOVEC LONG 
Guecenenseenyequeenenetend 
! ; Performance indicators ! 
Phase Page faults CPU Time Elapsed Time 
Initialization 31 00 300: 00.08 00:00:00.94 
zonnane processing 109 00:00:00.64 ee are 
6¢ BOD aR Bt 00:00:03.49 
aE table sort 0:00:00. ! Bo +20 :80-0 
Pass 2 61 0:00:00.6 80 :80:8 “7 
Symbol table output ¢ 0:00:00.03 B: 0:00.0 
Psect synopsis output 0:00:00.02 0:00:00.02 
Cross-reference output 8 00:00: 3-99 90:00:00.00 
Assembler run totals 29 00:00:02.26 00:00:11.12 


Tog work ing et Limit was 900 pages. 
byies (8 pages) of virtua memory were used to buffer the 
ok were 10 pages of symbol table space atieceiee to hold 17 non-local and 0 local symbols. 
330 source Lines were read in Pass 1, producing 11 object records in Pass 2. 
1 page of virtual memory was used to “define 1 macro. 
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T QRT 
He Macro Run Statistics 
0 GETS were required to define 0 macros. 
There were no errors, warnings or information messages. 


MACRO/ENABLE=SUPPRESSION/DI SABLE=(GLOBAL , TRACEBACK) /LIS=LIS$:MTHGSQRT/OBJ=OBJ$:MTHGSQRT MSRC$:MTHJACKET/UPDATE=(ENHS$:MTHJACKET) #MSRC 
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